import numpy as np
import sys
import os

gauss_line_weight = {1:np.array([2]),
                     2:np.array([1, 1]),
                     3:np.array([8/9, 5/9, 5/9]),
                     4:np.array([0.6521451548625462, 0.6521451548625462, 0.34785484513745385, 0.34785484513745385]),
                     5:np.array([0.5688888888888889, 0.47862867049936647, 0.47862867049936647, 0.23692688505618908, 0.23692688505618908])
                    }
gauss_line_point = {1:np.array([0.0]),
                    2:np.array([0.5773502691896258, -0.5773502691896258]),
                    3:np.array([0.0, 0.7745966692414834, -0.7745966692414834]),
                    4:np.array([0.3399810435848563, -0.3399810435848563, 0.8611363115940526, -0.8611363115940526]),
                    5:np.array([0.0, 0.5384693101056831, -0.5384693101056831, 0.906179845938664, -0.906179845938664])
                    }

gauss_tri_weight = {1:np.array([1]),
                    2:np.array([0.333333333333333, 0.333333333333333, 0.333333333333333]),
                    2.5:np.array([0.333333333333333, 0.333333333333333, 0.333333333333333]),
                    3:np.array([-0.562500000000000, 0.520833333333333, 0.520833333333333, 0.520833333333333]),
                    }

gauss_tri_point = {1:np.array([[1/3,1/3]]),
                   2:np.array([[0.166666666666667, 0.166666666666667],
                               [0.666666666666667, 0.166666666666667],
                               [0.166666666666667, 0.666666666666667]]),
                   2.5:np.array([[0.500000000000000, 0.000000000000000],
                               [0.000000000000000, 0.500000000000000],
                               [0.500000000000000, 0.500000000000000]]),
                   # 3:np.array([[0.333333333333333, 0.333333333333333],
                   #             [0.600000000000000, 0.200000000000000],
                   #             [0.600000000000000, 0.200000000000000],
                   #             [0.600000000000000, 0.200000000000000]]),
                   }


def gauss_line_int(f, x1, y1, x2, y2, pt=3):
    w = gauss_line_weight[pt]
    x = gauss_line_point[pt]
    
    L = np.sqrt((x1-x2)**2+(y1-y2)**2)

    ret = 0.5*L*np.sum(w*f(0.5*(x1-x2)*x+0.5*(x1+x2), 0.5*(y1-y2)*x+0.5*(y1+y2)))
    return ret


def gauss_tri_int(f, x1, y1, x2, y2, x3, y3, degree=2):
    
    area = 1/2*np.linalg.det(np.array([[x1,y1,1],
                                       [x2,y2,1],
                                       [x3,y3,1]]))
    
    x = lambda xi,eta: (x1-x3)*xi + (x2-x3)*eta + x3
    y = lambda xi,eta: (y1-y3)*xi + (y2-y3)*eta + y3
    
    w = gauss_tri_weight[degree]
    xi = gauss_tri_point[degree][:,0]
    eta = gauss_tri_point[degree][:,1]
    
    ret = area * w @ f(x(xi,eta),y(xi,eta))
    
    return ret

def show_cross_element(model):
    import numpy as np
    from mayavi import mlab

    x = model.Node[:,0]
    y = model.Node[:,1]
    z = np.zeros_like(x)
    triangles1 = model.Element
    triangles2 = model.Element[list(model.crossElement.keys()),:]

    mlab.triangular_mesh(x, y, z, triangles1, color=(1,1,1),representation='wireframe' )
    mlab.triangular_mesh(x, y, z, triangles2, color=(1,0,0),representation='surface' )

    for cluster in model.NeuCluster:
        i,_ = cluster[0]
        x = model.Node[i][0]
        y = model.Node[i][1]
        z = 0    
        mlab.points3d(x,y,z,1, color=(0,1,0),scale_factor=0.01)
        for i, j in cluster:
            x = model.Node[i][0]
            y = model.Node[i][1]
            z = 0
            u = model.Node[j][0] - x
            v = model.Node[j][1] - y
            w = 0
            mlab.quiver3d(x, y, z, u, v, w, scale_factor=1)

    for edge in model.crossEdge:
        i,j = edge
        x = [model.Node[i][0], model.Node[j][0]]
        y = [model.Node[i][1], model.Node[j][1]]
        z = [0,0]
        mlab.plot3d(x, y, z, color=(0,0,1),tube_radius=0.005)
            
            
    mlab.view(0,0)
    # print(mlab.view())
    mlab.show()

class MutePrints:
    def __enter__(self):
        self._original_stdout = sys.stdout
        sys.stdout = open(os.devnull, 'w')

    def __exit__(self, exc_type, exc_val, exc_tb):
        sys.stdout.close()
        sys.stdout = self._original_stdout


class ConvergenceRate:
    """docstring for ConvergenceRate"""
    def __init__(self, model, geometry_func):
        self.model = model
        self.geometry_func = geometry_func
        self.hs_efem = None
        self.hs_fem = None
        self.dofs_efem = None
        self.dofs_fem = None
        self.errs_efem = None
        self.errs_fem = None
        self.conds_fem = None
        self.conds_efem = None
        self.dof_slope_fem = 0
        self.size_slope_fem = 0
        self.dof_slope_efem = 0
        self.size_slope_efem =  0
        self.dof_slope_cre = 0
        self.size_slope_cre = 0
        self.query_efem = None
        self.query_fem = None
        self.rs = None
        self.rs_fem = None


    def convergence_efem(self, exact_ux, exact_uy, start=2, end=-1, num=6, display=False):
        self.hs_efem = np.logspace(start, end, num)
        self.errs_efem = []
        for h in self.hs_efem:
            self.model.Geometry = self.geometry_func(h)
            self.model.mesh()
            self.model.efem()
            self.errs_efem.append(self.model.error_efem(exact_ux, exact_uy))
        if display:
            self.plot_converage_rate(self.hs_efem, self.errs_efem)


    def convergence_fem(self, exact_ux, exact_uy, start=2, end=-1, num=6, display=False):
        self.hs_fem = np.logspace(start, end, num)
        self.errs_fem = []
        for h in self.hs_fem:
            self.model.Geometry = self.geometry_func(h)
            self.model.mesh()
            self.model.fem()
            self.errs_fem.append(self.model.error_fem(exact_ux, exact_uy))
        if display:
            self.plot_converage_rate(self.hs_fem, self.errs_fem)

    def plot_converage_rate(self, h, e):
        import matplotlib.pyplot as plt
        plt.figure()
        h = [1/x for x in h]
        plt.plot(h, e)
        plt.scatter(h, e, marker='o')
        plt.xscale("log")
        plt.yscale("log")
        plt.grid()
        plt.show()

    def convergence_all(self, exact_ux, exact_uy, start=2, end=-1, num=6, display=False, condition=False, query=False):
        self.hs_list = np.logspace(start, end, num)
        self.hs_fem = []
        self.dofs_fem = []
        self.errs_fem = []
        # self.hs_efem = np.logspace(start, end, num)
        self.hs_efem = []
        self.dofs_efem = []
        self.errs_efem = []
        self.cres = []
        self.dof_slope_fem = 0
        self.size_slope_fem = 0
        self.dof_slope_efem = 0
        self.size_slope_efem =  0
        self.dof_slope_cre = 0
        self.size_slope_cre = 0
        self.conds_fem = []
        self.conds_efem = []

        if self.model.numHole != 0:
            self.rs = []
            self.rs_fem = []

        if query:
            self.query_efem = []
            self.query_fem = []

        for h in self.hs_list:
            self.model.Geometry = self.geometry_func(h)
            self.model.mesh(display=display)
            self.model.fem()
            self.model.efem()
            self.dofs_fem.append(self.model.numDOF)
            self.dofs_efem.append(self.model.numDOF)
            self.hs_fem.append(1/np.sqrt(self.model.numDOF))
            self.hs_efem.append(1/np.sqrt(self.model.numDOF))
            self.errs_fem.append(self.model.error_fem(exact_ux, exact_uy))
            self.errs_efem.append(self.model.error_efem(exact_ux, exact_uy))
            self.cres.append(self.model.calculate_cre())
            if condition:
                self.conds_fem.append(np.linalg.cond(self.model.K))
                self.conds_efem.append(np.linalg.cond(self.model.F))
            if self.model.numHole != 0:
                self.rs_fem.append(self.model.rs_fem)
                self.rs.append([self.model.s[self.model.numNode:self.model.numNode+self.model.numHole]])

            if query:
                query_loc = np.where((self.model.Node[:,0]==query[0])*(self.model.Node[:,1]==query[1]))[0][0]
                # print(query_loc, self.model.Node[query_loc])
                if query_loc:
                    # neighbor_fem = []
                    # neighbor_efem = []

                    # for holeid, nodes in enumerate(self.model.Element):
                    #     if query_loc in nodes:
                    #         neighbor_efem.append(self.model.qh[holeid])
                    #         neighbor_fem.append(self.model.qh_fem[holeid])

                    # self.query_efem.append(np.mean(np.array(neighbor_efem),axis=0))
                    # self.query_fem.append(np.mean(np.array(neighbor_fem),axis=0))

                    for elemid, nodes in enumerate(self.model.Element):
                        if query_loc in nodes:
                            order = set(nodes)
                            order.remove(query_loc)
                            i = order.pop()
                            j = order.pop()
                            if np.allclose(self.model.Node[query_loc,0],self.model.Node[i,0]) or np.allclose(self.model.Node[query_loc,0],self.model.Node[j,0]):
                                self.query_efem.append(self.model.qh[elemid])
                                self.query_fem.append(self.model.qh_fem[elemid])
                                # print(elemid, self.model.Element[elemid])
                                # print(self.model.Node[self.model.Element[elemid]])
                                break
                else:
                    print("query point not found!")

            print('Calculation complete for mesh size:{:.4f}'.format(h))
        
        self.query_efem = np.array(self.query_efem)
        self.query_fem = np.array(self.query_fem)

        A = np.c_[np.log(self.dofs_fem), np.ones_like(self.dofs_fem)]
        y = np.log(self.errs_fem)
        k, b = np.linalg.inv(A.T@A) @ (A.T@y)
        self.dof_slope_fem = k

        A = np.c_[np.log(self.hs_fem), np.ones_like(self.hs_fem)]
        y = np.log(self.errs_fem)
        k, b = np.linalg.inv(A.T@A) @ (A.T@y)
        self.size_slope_fem = k

        A = np.c_[np.log(self.dofs_efem), np.ones_like(self.dofs_efem)]
        y = np.log(self.errs_efem)
        k, b = np.linalg.inv(A.T@A) @ (A.T@y)
        self.dof_slope_efem = k

        A = np.c_[np.log(self.hs_efem), np.ones_like(self.hs_efem)]
        y = np.log(self.errs_efem)
        k, b = np.linalg.inv(A.T@A) @ (A.T@y)
        self.size_slope_efem =  k     

        A = np.c_[np.log(self.dofs_efem), np.ones_like(self.dofs_efem)]
        y = np.log(self.cres)
        k, b = np.linalg.inv(A.T@A) @ (A.T@y)
        self.dof_slope_cre = k

        A = np.c_[np.log(self.hs_efem), np.ones_like(self.hs_efem)]
        y = np.log(self.cres)
        k, b = np.linalg.inv(A.T@A) @ (A.T@y)
        self.size_slope_cre =  k   


        # if condition:
        #     print("Convergence rate calculation is complete!")
        #     print("Condition number:")
        #     print("Size \t|  FEM \t| EFEM  ")
        #     for h, fem, efem in zip(self.hs_fem, self.conds_fem, self.conds_efem):
        #         print('{:.6f} \t| {:.6e} \t| {:.6e} '.format(h, fem, efem))

    def convergence_all_without_exact(self, start=2, end=-1, num=6, display=False):
        self.hs_list = np.logspace(start, end, num)
        self.hs_fem = []
        self.dofs_fem = []
        # self.hs_efem = np.logspace(start, end, num)
        self.hs_efem = []
        self.dofs_efem = []
        self.cres = []
        self.dof_slope_cre = 0
        self.size_slope_cre = 0

        if self.model.numHole != 0:
            self.rs = []
            self.rs_fem = []
        # if condition:
        #     self.conds_fem = []
        #     self.conds_efem = []
        # if query:
        #     self.query_efem = []
        #     self.query_fem = []


        for h in self.hs_list:
            self.model.Geometry = self.geometry_func(h)
            self.model.mesh(display=display)
            self.model.fem()
            self.model.efem()
            self.dofs_fem.append(self.model.numDOF)
            self.dofs_efem.append(self.model.numDOF)
            self.hs_fem.append(1/np.sqrt(self.model.numDOF))
            self.hs_efem.append(1/np.sqrt(self.model.numDOF))
            self.cres.append(self.model.calculate_cre())
            print('Calculation complete for mesh size:{:.4f}'.format(h))
            if self.model.numHole != 0:
                self.rs_fem.append(self.model.rs_fem)
                self.rs.append([self.model.s[self.model.numNode:self.model.numNode+self.model.numHole]])


        A = np.c_[np.log(self.dofs_efem), np.ones_like(self.dofs_efem)]
        y = np.log(self.cres)
        k, b = np.linalg.inv(A.T@A) @ (A.T@y)
        self.dof_slope_cre = k

        A = np.c_[np.log(self.hs_efem), np.ones_like(self.hs_efem)]
        y = np.log(self.cres)
        k, b = np.linalg.inv(A.T@A) @ (A.T@y)
        self.size_slope_cre =  k 


    def plot_multiple_convergence_rate(self, which='size', file_name=None, font_size=18, fem_kwargs={}, efem_kwargs={}, fem_aux_kwargs={}, efem_aux_kwargs={}):
        if (self.hs_efem is None) or (self.hs_fem is None):
            print("Please call convergence_all() method first!")
            return

        hs_group = [self.hs_fem, self.hs_efem]
        dofs_group = [self.dofs_fem, self.dofs_efem]
        errs_group = [self.errs_fem, self.errs_efem]
        dofs_slope_group = [self.dof_slope_fem, self.dof_slope_efem]
        size_slope_group = [self.size_slope_fem, self.size_slope_efem]

        kwargs_group = [fem_kwargs, efem_kwargs] 
        aux_kwargs_group = [fem_aux_kwargs, efem_aux_kwargs] 


        import matplotlib.pyplot as plt
        import matplotlib

        matplotlib.rc('text', usetex=True)
        matplotlib.rc('text.latex', preamble=r'\usepackage{amsmath}')

        fig, ax = plt.subplots(figsize=(4,6), dpi=300)

        if which == "size":
            for h, e, s, kwargs in zip(hs_group, errs_group, size_slope_group, kwargs_group):
                kwargs['label'] = kwargs['label'] + '(Slope: {:.2f})'.format(s)
                ax.plot(h, e, **kwargs)

            # for h, aux_kwargs in zip(hs_group, aux_kwargs_group):
            #     aux = [aux_kwargs["C"]*x**-1 for x in h]
            #     del aux_kwargs["C"]
            #     ax.plot(h, aux, **aux_kwargs)
            # ax.set_xlabel(r"Inverse of mesh size, $\frac{1}{h}$")
        elif which == "dofs":
            for d, e, s, kwargs in zip(dofs_group, errs_group, dofs_slope_group, kwargs_group):
                kwargs['label'] = kwargs['label'].format(s)
                ax.plot(d, e, **kwargs)

            # for d, aux_kwargs in zip(dofs_group, aux_kwargs_group):
            #     aux = [aux_kwargs["C"]*x**-1 for x in d]
            #     del aux_kwargs["C"]
            #     ax.plot(d, aux, **aux_kwargs)
            # ax.set_xlabel(r"Degrees of freedom, $n$")
        else:
            print("The parameter 'which' should be 'size' or 'dofs'!")
            return    

        ax.tick_params(axis="both", which="both", direction="in")    
        ax.tick_params(axis='both', which='major', labelsize=font_size)
        ax.tick_params(axis='both', which='major', length=8)
        ax.tick_params(axis='both', which='major', width=1)
        ax.tick_params(axis='both', which='minor', length=4)
        ax.tick_params(axis='both', which='minor', width=1)

        # ax.tick_params(axis = 'both', which = 'minor', labelsize = 24)

        ax.set_xscale("log")
        ax.set_yscale("log")
        # ax.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        
        # ax.set_ylabel(r"Error in norm, $e$")
        # ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        # ax.xaxis.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%d'))
        ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.LogFormatter(labelOnlyBase=True, minor_thresholds=(0.5,0.2)))

        # ax.grid()
        plt.setp(ax.spines.values(), linewidth=2)
        plt.legend(fontsize='x-large')
        if file_name is None:
            plt.savefig("multiple_convergence.png")
        else:
            plt.savefig(file_name)
        plt.show()
        

class CompareResult:
    """docstring for CompareResult"""
    def __init__(self, model=None, geometry_func=None):
        self.model = model
        self.geometry_func = geometry_func
        self.triangles = None
        self.x = None
        self.y = None

    # Final version.
    def cmp_table(self, cr_model, file_name='table.pkl'):
        import pickle

        kv = {'Size'         : (cr_model.hs_efem, '{:.4f}'),
              'DOFs'         : (cr_model.dofs_efem, '{:d}'),
              '||qh-duh||_q' : (cr_model.cres, '{:.4f}'),
              }

        if cr_model.errs_efem:
            kv['||q-qh||_q'] = (cr_model.errs_efem, '{:.4f}')
            kv['||u-uh||_u'] = (cr_model.errs_fem, '{:.4f}')

        if cr_model.conds_efem:
            kv['cond(F)'] = (cr_model.conds_efem, '{:.2e}')
            kv['cond(K)'] = (cr_model.conds_fem, '{:.2e}')

        # if cr_model.query_efem.any():
        #     kv['qx(efem)'] = (cr_model.query_efem[:,0], '{:.4f}')
        #     kv['qx(fem)']  = ( cr_model.query_fem[:,0], '{:.4f}')
        #     kv['qy(efem)'] = (cr_model.query_efem[:,1], '{:.4f}')
        #     kv['qy(fem)']  = ( cr_model.query_fem[:,1], '{:.4f}')

        if cr_model.rs:
            rs = np.array([x[0] for x in cr_model.rs])
            rs_fem = np.array(cr_model.rs_fem)
            for holeid in range(rs.shape[1]):
                kv['r{}(efem)'.format(holeid+1)] = (rs[:,holeid], '{:.4f}')
                kv['r{}(fem)'.format(holeid+1)]  = (rs_fem[:,holeid], '{:.4f}')         

        order_key = [key for key in kv]
        heading = "\t".join([ key for key in kv])
        # fmt = "\t".join([ kv[key][1] for key in kv])
        fmt = "\t".join([ '{}' for key in kv])
        groups = [kv[key][0] for key in kv]
        # print(groups)
        print(heading)
        
        for item in zip(*groups):
            print(fmt.format(*item))

        with open(file_name,"wb") as f:
            pickle.dump(kv,f)



        # if cr_model.conds_fem:
        #     print("   Size   \t|  cond_F \t| cond_K \t| ||q-qh||_q\t | ||u-uh||_u\t | ||qh-duh||_q \t|")
        #     for h, cond_F, cond_K, err_efem, err_fem, cre in zip(cr_model.hs_fem, cr_model.conds_efem, cr_model.conds_fem ,
        #                                  cr_model.errs_efem, cr_model.errs_fem, cr_model.cres):
        #         print('{:.6f}\t|{:.6e}\t|{:.6e}\t|{:.6f}\t|{:.6f}\t|{:.6f}\t|'.format(h, cond_K, cond_F, err_efem, err_fem, cre))        
        # else:
        #     print("   Size  \t| ||q-qh||_q\t| ||u-uh||_u\t| ||qh-duh||_q\t|")
        #     print("-------------------------------------------------------------------------")
        #     for h, err_efem, err_fem, cre in zip(cr_model.hs_fem, cr_model.errs_efem, cr_model.errs_fem, cr_model.cres):
        #         print('{:.4f}\t|{:.4f}\t|{:.4f}\t|{:.4f}\t|'.format(h, err_efem, err_fem, cre))        

    # Final version.
    def cmp_point_convergence_rate(self, cr_model, ux, uy, pt, which='size', direction='x', file_name="multiple_convergence.png", font_size=18, fem_kwargs={}, efem_kwargs={}):
        if (cr_model.hs_efem is None) or (cr_model.hs_fem is None):
            print("Please call convergence_all() method first!")
            return

        import matplotlib.pyplot as plt
        import matplotlib

        # matplotlib.rc('text', usetex=True)
        # matplotlib.rc('text.latex', preamble=r'\usepackage{amsmath}')

        fig, ax = plt.subplots(dpi=300)
        # ax.set_aspect('equal', 'box')
        # ax.set_aspect('auto')
        ax.set_box_aspect(1)

        if direction == "x":
            true_value = ux(*pt)
            query_efem = cr_model.query_efem[:,0]
            query_fem = cr_model.query_fem[:,0]
        elif direction == "y":
            true_value = uy(*pt)
            query_efem = cr_model.query_efem[:,1]
            query_fem = cr_model.query_fem[:,1]
        else:
            print("direction should be 'x' or 'y'!")
            return

        if which == "size":
            A = np.c_[np.log(cr_model.hs_efem[2::]), np.ones_like(cr_model.hs_efem[2::])]
            y1 = np.log(np.abs((query_efem[2::]-true_value)/true_value))
            y2 = np.log(np.abs((query_fem[2::]-true_value)/true_value))
            slope1, _ = np.linalg.inv(A.T@A) @ (A.T@y1)   
            slope2, _ = np.linalg.inv(A.T@A) @ (A.T@y2)   

            efem_kwargs['label'] = efem_kwargs['label'] + '(Slope: {:.2f})'.format(slope1)
            fem_kwargs['label'] = fem_kwargs['label'] + '(Slope: {:.2f})'.format(slope2)
            ax.plot(cr_model.hs_efem, np.abs((query_efem-true_value)/true_value), **efem_kwargs)
            ax.plot(cr_model.hs_fem, np.abs((query_fem-true_value)/true_value),**fem_kwargs)
            # ax.plot(cr_model.hs_efem, np.abs(query_efem-true_value), **efem_kwargs)
            # ax.plot(cr_model.hs_fem, np.abs(query_fem-true_value),**fem_kwargs)
        elif which == "dofs":
            A = np.c_[np.log(cr_model.dofs_efem), np.ones_like(cr_model.dofs_efem)]
            y1 = np.log(np.abs((query_efem-true_value)/true_value))
            y2 = np.log(np.abs((query_fem-true_value)/true_value))
            slope1, _ = np.linalg.inv(A.T@A) @ (A.T@y1)   
            slope2, _ = np.linalg.inv(A.T@A) @ (A.T@y2)   

            efem_kwargs['label'] = efem_kwargs['label'] + '(Slope: {:.2f})'.format(slope1)
            fem_kwargs['label'] = fem_kwargs['label'] + '(Slope: {:.2f})'.format(slope2)
            ax.plot(cr_model.dofs_efem, np.abs((query_efem-true_value)/true_value), **efem_kwargs)
            ax.plot(cr_model.dofs_fem, np.abs((query_fem-true_value)/true_value), **fem_kwargs)
        else:
            print("The parameter 'which' should be 'size' or 'dofs'!")
            return    

        ax.tick_params(axis="both", which="both", direction="in")    
        ax.tick_params(axis='both', which='major', labelsize=font_size)
        ax.tick_params(axis='both', which='major', length=8)
        ax.tick_params(axis='both', which='major', width=1)
        ax.tick_params(axis='both', which='minor', length=4)
        ax.tick_params(axis='both', which='minor', width=1)
        # ax.tick_params(axis = 'both', which = 'minor', labelsize = 24)
        ax.set_xscale("log")
        ax.set_yscale("log")
        # ax.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        # ax.set_ylabel(r"Error in norm, $e$")
        # ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        # ax.xaxis.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%d'))
        ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.LogFormatter(labelOnlyBase=True, minor_thresholds=(0.5,0.2)))

        # ax.set_xlim([0.0085,0.04])
        # ax.set_ylim([0.05,20])
        plt.setp(ax.spines.values(), linewidth=2)
        # plt.legend(fontsize='large')
        plt.legend(fontsize='large', frameon=False)
        plt.tight_layout()

        plt.savefig(file_name)
        plt.show()        

    def cmp_resultant_force(self, cr_model, q0n, holeid=0, file_name="hole_resultant_force.png", font_size=18, fem_kwargs={}, efem_kwargs={}):
        import matplotlib.pyplot as plt
        import matplotlib

        matplotlib.rc('text', usetex=True)
        matplotlib.rc('text.latex', preamble=r'\usepackage{amsmath}')

        fig, ax = plt.subplots(dpi=300)
        # ax.set_aspect('equal', 'box')
        # ax.set_aspect('auto')
        ax.set_box_aspect(1)

        rs = np.array([x[0] for x in cr_model.rs])
        rs_fem = np.array(cr_model.rs_fem)
        h = 1/np.sqrt(cr_model.dofs_fem)

        ax.plot(h, rs[:,holeid]+q0n[holeid], **efem_kwargs)
        ax.plot(h, rs_fem[:,holeid], **fem_kwargs)

        ax.tick_params(axis="both", which="both", direction="in")    
        ax.tick_params(axis='both', which='major', labelsize=font_size)
        ax.tick_params(axis='both', which='major', length=8)
        ax.tick_params(axis='both', which='major', width=1)
        ax.tick_params(axis='both', which='minor', length=4)
        ax.tick_params(axis='both', which='minor', width=1)
        # ax.tick_params(axis = 'both', which = 'minor', labelsize = 24)
        ax.set_xscale("log")
        # ax.set_yscale("log")
        # ax.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        # ax.set_ylabel(r"Error in norm, $e$")
        # ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        # ax.xaxis.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%d'))
        ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.LogFormatter(labelOnlyBase=True, minor_thresholds=(0.5,0.2)))
        ax.set_xlim([0.008,0.105])
        # ax.grid()
        # ax.set_xlim([0.005,0.5])
        # ax.set_ylim([0.05,20])
        plt.setp(ax.spines.values(), linewidth=2)
        # plt.legend(fontsize='large')
        plt.legend(fontsize='large', frameon=False)
        plt.tight_layout()

        plt.savefig(file_name)
        plt.show()    

    def cmp_matrix(self, mesh_size, display=False, file_name='cmp.png',cmap='viridis'):
        self.model.Geometry = self.geometry_func(mesh_size)
        self.model.mesh(display=display)
        self.model.fem()
        self.model.efem()

        import matplotlib.pyplot as plt
        from matplotlib import cm
        import matplotlib.patches as mpatches

        fig = plt.figure(dpi=300)
        axs = fig.subplots(1,2)

        matrix = [self.model.KU, self.model.FU]
        for i, ax in enumerate(axs):
            ax.matshow((matrix[i] != 0) * 1,cmap=cmap)
            ax.axis('off')

        cf = cm.get_cmap(cmap,2)
        colors = {1:cf(0), 2:cf(1)}
        labels = {1:'Zero', 2:'Non-zero'}
        patches =[mpatches.Patch(color=colors[i],label=labels[i]) for i in colors]

        plt.legend(handles=patches, bbox_to_anchor=(0,1.15), loc="upper left", mode="expand", ncol=2,borderaxespad=0, frameon=False)
        plt.savefig(file_name)
        plt.show()
        print(np.allclose(self.model.KU, self.model.FU))

    def cmp_field(self, ux, uy, mesh_size, display=False, file_name='field.svg', **kwargs):
        self.model.Geometry = self.geometry_func(mesh_size)
        self.model.mesh(display=display)
        # self.model.fem(display=display)
        self.model.efem(display=display)

        self.triangles = self.model.Element.copy()
        self.x = self.model.Node[:,0].copy()
        self.y = self.model.Node[:,1].copy()
        self.ux = ux(self.x, self.y)
        self.uy = uy(self.x, self.y)

        mm = np.concatenate([self.ux.flatten(),self.uy.flatten(),self.model.qh[:,0].flatten(),self.model.qh[:,1]])
        vmin = np.min(mm)
        vmax = np.max(mm)

        import matplotlib.pyplot as plt
        fig = plt.figure(dpi=300)
        axs = fig.subplots(2,2, sharex=True, sharey=True)

        tc = axs[0,0].tricontourf(self.x, self.y, self.triangles, self.ux, levels=256, vmin=vmin, vmax=vmax, **kwargs)
        axs[0,0].set_aspect('equal','box')
        axs[0,0].tick_params(bottom=False)
        # axs[0,0].tick_params(left=False, bottom=False)

        axs[1,0].tripcolor(self.x, self.y, self.triangles, facecolors=self.model.qh[:,0], vmin=vmin, vmax=vmax, **kwargs)
        axs[1,0].set_aspect('equal','box')

        axs[0,1].tricontourf(self.x, self.y, self.triangles, self.uy, levels=256, vmin=vmin, vmax=vmax, **kwargs)
        axs[0,1].set_aspect('equal','box')
        axs[0,1].tick_params(left=False, bottom=False)

        axs[1,1].tripcolor(self.x, self.y, self.triangles, facecolors=self.model.qh[:,1], vmin=vmin, vmax=vmax, **kwargs)
        axs[1,1].set_aspect('equal','box')
        axs[1,1].tick_params(left=False)


        axs[0,0].set_yticks([0,0.5,1.0])
        axs[1,0].set_yticks([0,0.5,1.0])

        fig.colorbar(tc, ax=axs.ravel().tolist())

        # fig.tight_layout()
        plt.savefig(file_name)
        plt.show()
        # self.model.uh
        # self.model.qh

    def cmp_two_method_results(self, mesh_size, display=False, file_name='two_method_results.svg', **kwargs):
        from copy import deepcopy
        self.model.Geometry = self.geometry_func(mesh_size)
        self.model.mesh(display=display)
        self.model.efem(display=display)
        self.model.fem(display=display)


        triangles = self.model.Element.copy()
        x = self.model.Node[:,0].copy()
        y = self.model.Node[:,1].copy()

        n = len(x)
        qh_node = np.zeros((n,2))
        qh_node_fem = np.zeros((n,2))

        for nodeid in range(n):
            elemids = np.where((self.model.Element==nodeid).any(axis=1)==True)
            qh_node[nodeid] = np.mean(self.model.qh[elemids],axis=0)
            qh_node_fem[nodeid] = np.mean(self.model.qh_fem[elemids],axis=0)

        mm1 = np.concatenate([qh_node[:,0].flatten(),qh_node_fem[:,0].flatten()])
        vmin1 = np.min(mm1)
        vmax1 = np.max(mm1)

        mm2 = np.concatenate([qh_node[:,1].flatten(),qh_node_fem[:,1].flatten()])
        vmin2 = np.min(mm2)
        vmax2 = np.max(mm2)

        # mm1 = np.concatenate([self.model.qh[:,0].flatten(),self.model.qh_fem[:,0].flatten()])
        # vmin1 = np.min(mm1)
        # vmax1 = np.max(mm1)

        # mm2 = np.concatenate([self.model.qh[:,1].flatten(),self.model.qh_fem[:,1].flatten()])
        # vmin2 = np.min(mm2)
        # vmax2 = np.max(mm2)

        import matplotlib.pyplot as plt

        def fmt(x):
            s = f"{x:.2f}"
            return rf"{s}" if plt.rcParams["text.usetex"] else f"{s}"

        levels = 15

        fig, axs = plt.subplots(nrows=2,ncols=2,dpi=300,sharex=True, sharey=True, figsize=(4,3))
        fig.subplots_adjust(wspace=0.1)
        # tc = axs[0,0].tripcolor(x, y, triangles, facecolors=self.model.qh[:,0], vmin=vmin1, vmax=vmax1, **kwargs)
        
        tc = axs[0,0].tricontour(x, y, triangles, qh_node[:,0], linewidths=0.2, colors='k', levels=levels, vmin=vmin1, vmax=vmax1)
        tf = axs[0,0].tricontourf(x, y, triangles, qh_node[:,0], levels=levels, vmin=vmin1, vmax=vmax1, **kwargs)
        # axs[0,0].clabel(tc, tc.levels, inline=True, fmt=fmt, fontsize=6, colors='white')

        axs[0,0].set_aspect('equal','box')
        axs[0,0].tick_params(bottom=False)
        axs[0,0].spines["top"].set_visible(False)
        axs[0,0].spines["right"].set_visible(False)
        axs[0,0].spines["bottom"].set_visible(False)
        # axs[0,0].set_frame_on(False)
        # axs[0,0].axis('off')

        # axs[0,1].tripcolor(x, y, triangles, facecolors=self.model.qh_fem[:,0], vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,1].tricontour(x, y, triangles, qh_node_fem[:,0], linewidths=0.2, colors='k', levels=levels, vmin=vmin1, vmax=vmax1)
        axs[0,1].tricontourf(x, y, triangles, qh_node_fem[:,0], levels=levels, vmin=vmin1, vmax=vmax1, **kwargs)

        axs[0,1].set_aspect('equal','box')
        axs[0,1].tick_params(left=False, bottom=False)
        # axs[0,1].set_frame_on(False)
        axs[0,1].axis('off')

        fig.colorbar(tf, ax=axs[0,:].flat)

        # tc = axs[1,0].tripcolor(x, y, triangles, facecolors=self.model.qh[:,1], vmin=vmin2, vmax=vmax2, **kwargs)

        tc = axs[1,0].tricontour(x, y, triangles, qh_node[:,1], linewidths=0.2, colors='k', levels=levels, vmin=vmin1, vmax=vmax1)
        tf = axs[1,0].tricontourf(x, y, triangles, qh_node[:,1], levels=levels, vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,0].set_aspect('equal','box')
        axs[1,0].spines["top"].set_visible(False)
        axs[1,0].spines["right"].set_visible(False)

        # axs[1,1].tripcolor(x, y, triangles, facecolors=self.model.qh_fem[:,1], vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,1].tricontour(x, y, triangles, qh_node_fem[:,1], linewidths=0.2, colors='k', levels=levels, vmin=vmin1, vmax=vmax1)
        axs[1,1].tricontourf(x, y, triangles, qh_node_fem[:,1], levels=levels, vmin=vmin2, vmax=vmax2, **kwargs)        

        axs[1,1].set_aspect('equal','box')
        axs[1,1].tick_params(left=False)
        axs[1,1].spines["top"].set_visible(False)
        axs[1,1].spines["right"].set_visible(False)
        axs[1,1].spines["left"].set_visible(False)

        fig.colorbar(tf, ax=axs[1,:].flat)

        axs[0,0].set_yticks([0,0.5,1.0])
        axs[1,0].set_yticks([0,0.5,1.0])
        axs[1,0].set_xticks([0,0.5,1.0])
        axs[1,1].set_xticks([0,0.5,1.0])

        axs[0,0].set_xlim([0,1])
        axs[0,0].set_ylim([0,1])
        axs[1,0].set_ylim([0,1])
        # for ax in axs:
        #     ax.set_xlim([0,1])
        #     ax.set_ylim([0,1])

        # fig.tight_layout()
        plt.savefig(file_name)
        plt.show()
        # self.model.uh
        # self.model.qh

    # 为检查审稿人提出的疑问所写的检查程序
    def _flux(self, ux, uy, q0_x, q0_y, mesh_size, direction=0, display=False, file_name='two_mesh_field.svg', **kwargs ):
        from copy import deepcopy
        self.model.Geometry = self.geometry_func(mesh_size)
        self.model.mesh(display=display)
        self.model.efem(display=display)
        q0 = [q0_x,q0_y]
        tq0 = q0[direction]
        qh = deepcopy(self.model.qh0)
        x = self.model.Node[:,0].copy()
        y = self.model.Node[:,1].copy()
        tri = self.model.Element.copy()
        print(tri.shape)
        vmin = np.min(qh[:,direction].flatten())
        vmax = np.max(qh[:,direction].flatten())

        import matplotlib.pyplot as plt
        fig, ax = plt.subplots(dpi=150,figsize=(6,3))
        for i, t in enumerate(tri):
            tx = x[t]
            ty = y[t]
            ax.tricontourf(tx, ty, [[0,1,2]], qh[i,direction]+tq0(tx,ty), vmin=vmin, vmax=vmax,levels=50, **kwargs)
            # ax.tripcolor(tx, ty, [[0,1,2]], qh[i,direction]+tq0(tx,ty), vmin=vmin, vmax=vmax, shading='flat', **kwargs)
        
        ax.set_aspect('equal','box')
        plt.show()
        return x,y,tri,qh
    # 为检查审稿人提出的疑问所写的检查程序
    def cmp_two_mesh_field_new(self, ux, uy, mesh_size1, mesh_size2, q0_x, q0_y, display=False, file_name='two_mesh_field.svg', **kwargs):
        from copy import deepcopy
        self.model.Geometry = self.geometry_func(mesh_size1)
        self.model.mesh(display=display)
        self.model.efem(display=display)
        first_qh = deepcopy(self.model.qh0)
        _first_qh = deepcopy(self.model.qh)
        first_x = self.model.Node[:,0].copy()
        first_y = self.model.Node[:,1].copy()
        first_tri = self.model.Element.copy()


        self.model.Geometry = self.geometry_func(mesh_size2)
        self.model.mesh(display=display)
        self.model.efem(display=display)
        second_qh = deepcopy(self.model.qh)
        second_x = self.model.Node[:,0].copy()
        second_y = self.model.Node[:,1].copy()
        second_tri = self.model.Element.copy()

        self.triangles = self.model.Element.copy()
        self.x = self.model.Node[:,0].copy()
        self.y = self.model.Node[:,1].copy()
        self.ux = ux(self.x, self.y)
        self.uy = uy(self.x, self.y)

        mm1 = np.concatenate([self.ux.flatten(),
                             _first_qh[:,0].flatten(),
                             second_qh[:,0].flatten(),])
        vmin1 = np.min(mm1)
        vmax1 = np.max(mm1)

        mm2 = np.concatenate([self.uy.flatten(),
                             _first_qh[:,1],
                             second_qh[:,1]])

        vmin2 = np.min(mm2)
        vmax2 = np.max(mm2)

        import matplotlib.pyplot as plt
        # fig = plt.figure(dpi=300)
        # axs = fig.subplots(2,3, sharex=True, sharey=True, constrained_layout=True)
        # fig.subplots_adjust(hspace=0.1)


        fig, axs = plt.subplots(nrows=2,ncols=3,dpi=300,sharex=True, sharey=True, figsize=(6,3))
        fig.subplots_adjust(wspace=0.1)
        tc1 = axs[0,0].tricontourf(self.x, self.y, self.triangles, self.ux, levels=256, vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,0].set_aspect('equal','box')
        axs[0,0].tick_params(bottom=False)
        # axs[0,0].tick_params(left=False, bottom=False)
        tc2 = axs[1,0].tricontourf(self.x, self.y, self.triangles, self.uy, levels=256, vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,0].set_aspect('equal','box')

        for i, t in enumerate(first_tri):
            tx = first_x[t]
            ty = first_y[t]
            axs[0,1].tricontourf(tx, ty, [[0,1,2]], first_qh[i,0]+q0_x(tx,ty), vmin=vmin1, vmax=vmax1,levels=100, **kwargs)
        # axs[0,1].tripcolor(first_x, first_y, first_tri, facecolors=first_qh[:,0], vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,1].set_aspect('equal','box')
        axs[0,1].tick_params(left=False, bottom=False)

        for i, t in enumerate(first_tri):
            tx = first_x[t]
            ty = first_y[t]
            axs[1,1].tricontourf(tx, ty, [[0,1,2]], first_qh[i,1]+q0_y(tx,ty), vmin=vmin2, vmax=vmax2,levels=100, **kwargs)
        # axs[1,1].tripcolor(first_x, first_y, first_tri, facecolors=first_qh[:,1], vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,1].set_aspect('equal','box')
        axs[1,1].tick_params(left=False)

        # for i, t in enumerate(second_tri):
        #     tx = second_x[t]
        #     ty = second_y[t]
        #     axs[0,2].tricontourf(tx, ty, [[0,1,2]], second_qh[i,0]+q0_x(tx,ty), vmin=vmin1, vmax=vmax1,levels=100, **kwargs)
        axs[0,2].tripcolor(second_x, second_y, second_tri, facecolors=second_qh[:,0], vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,2].set_aspect('equal','box')
        axs[0,2].tick_params(left=False, bottom=False)

        # for i, t in enumerate(second_tri):
        #     tx = second_x[t]
        #     ty = second_y[t]
        #     axs[1,2].tricontourf(tx, ty, [[0,1,2]], second_qh[i,1]+q0_y(tx,ty), vmin=vmin2, vmax=vmax2,levels=100, **kwargs)
        axs[1,2].tripcolor(second_x, second_y, second_tri, facecolors=second_qh[:,1], vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,2].set_aspect('equal','box')
        axs[1,2].tick_params(left=False)

        axs[0,0].set_yticks([0,0.5,1.0])
        axs[1,0].set_yticks([0,0.5,1.0])
        axs[1,0].set_xticks([0,0.5,1.0])
        axs[1,1].set_xticks([0,0.5,1.0])
        axs[1,2].set_xticks([0,0.5,1.0])


        fig.colorbar(tc1, ax=axs[0,:].flat)

        fig.colorbar(tc2, ax=axs[1,:].flat)

        # fig.tight_layout()
        plt.savefig(file_name)
        plt.show()
        # self.model.uh
        # self.model.qh        


    def _cmp_two_mesh_field(self, ux, uy, mesh_size1, mesh_size2, display=False, file_name='two_mesh_field.svg', **kwargs):
        from copy import deepcopy
        self.model.Geometry = self.geometry_func(mesh_size1)
        self.model.mesh(display=display)
        self.model.efem(display=display)
        first_qh = deepcopy(self.model.qh)
        first_x = self.model.Node[:,0].copy()
        first_y = self.model.Node[:,1].copy()
        first_tri = self.model.Element.copy()


        self.model.Geometry = self.geometry_func(mesh_size2)
        self.model.mesh(display=display)
        self.model.efem(display=display)
        second_qh = deepcopy(self.model.qh)
        second_x = self.model.Node[:,0].copy()
        second_y = self.model.Node[:,1].copy()
        second_tri = self.model.Element.copy()

        self.triangles = self.model.Element.copy()
        self.x = self.model.Node[:,0].copy()
        self.y = self.model.Node[:,1].copy()
        self.ux = ux(self.x, self.y)
        self.uy = uy(self.x, self.y)

        mm1 = np.concatenate([self.ux.flatten(),
                             first_qh[:,0].flatten(),
                             second_qh[:,0].flatten(),])
        vmin1 = np.min(mm1)
        vmax1 = np.max(mm1)

        mm2 = np.concatenate([self.uy.flatten(),
                             first_qh[:,1],
                             second_qh[:,1]])

        vmin2 = np.min(mm2)
        vmax2 = np.max(mm2)

        import matplotlib.pyplot as plt
        # fig = plt.figure(dpi=300)
        # axs = fig.subplots(2,3, sharex=True, sharey=True, constrained_layout=True)
        # fig.subplots_adjust(hspace=0.1)


        fig, axs = plt.subplots(nrows=2,ncols=3,dpi=300,sharex=True, sharey=True, figsize=(6,3))
        fig.subplots_adjust(wspace=0.1)
        tc1 = axs[0,0].tricontourf(self.x, self.y, self.triangles, self.ux, levels=256, vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,0].set_aspect('equal','box')
        axs[0,0].tick_params(bottom=False)
        # axs[0,0].tick_params(left=False, bottom=False)
        tc2 = axs[1,0].tricontourf(self.x, self.y, self.triangles, self.uy, levels=256, vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,0].set_aspect('equal','box')

        axs[0,1].tripcolor(first_x, first_y, first_tri, facecolors=first_qh[:,0], vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,1].set_aspect('equal','box')
        axs[0,1].tick_params(left=False, bottom=False)

        axs[1,1].tripcolor(first_x, first_y, first_tri, facecolors=first_qh[:,1], vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,1].set_aspect('equal','box')
        axs[1,1].tick_params(left=False)

        axs[0,2].tripcolor(second_x, second_y, second_tri, facecolors=second_qh[:,0], vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,2].set_aspect('equal','box')
        axs[0,2].tick_params(left=False, bottom=False)

        axs[1,2].tripcolor(second_x, second_y, second_tri, facecolors=second_qh[:,1], vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,2].set_aspect('equal','box')
        axs[1,2].tick_params(left=False)

        axs[0,0].set_yticks([0,0.5,1.0])
        axs[1,0].set_yticks([0,0.5,1.0])
        axs[1,0].set_xticks([0,0.5,1.0])
        axs[1,1].set_xticks([0,0.5,1.0])
        axs[1,2].set_xticks([0,0.5,1.0])


        fig.colorbar(tc1, ax=axs[0,:].flat)

        fig.colorbar(tc2, ax=axs[1,:].flat)

        # fig.tight_layout()
        plt.savefig(file_name)
        plt.show()
        # self.model.uh
        # self.model.qh


    def cmp_two_mesh_field(self, ux, uy, mesh_size1, mesh_size2, display=False, file_name='two_mesh_field.svg', **kwargs):
        from copy import deepcopy
        self.model.Geometry = self.geometry_func(mesh_size1)
        self.model.mesh(display=display)
        self.model.efem(display=display)
        first_qh = deepcopy(self.model.qh)
        first_x = self.model.Node[:,0].copy()
        first_y = self.model.Node[:,1].copy()
        first_tri = self.model.Element.copy()

        n = len(first_x)
        first_qh_node = np.zeros((n,2))

        for nodeid in range(n):
            elemids = np.where((self.model.Element==nodeid).any(axis=1)==True)
            first_qh_node[nodeid] = np.mean(self.model.qh[elemids],axis=0)



        self.model.Geometry = self.geometry_func(mesh_size2)
        self.model.mesh(display=display)
        self.model.efem(display=display)
        second_qh = deepcopy(self.model.qh)
        second_x = self.model.Node[:,0].copy()
        second_y = self.model.Node[:,1].copy()
        second_tri = self.model.Element.copy()

        n = len(second_x)
        second_qh_node = np.zeros((n,2))

        for nodeid in range(n):
            elemids = np.where((self.model.Element==nodeid).any(axis=1)==True)
            second_qh_node[nodeid] = np.mean(self.model.qh[elemids],axis=0)


        self.triangles = self.model.Element.copy()
        self.x = self.model.Node[:,0].copy()
        self.y = self.model.Node[:,1].copy()
        self.ux = ux(self.x, self.y)
        self.uy = uy(self.x, self.y)

        mm1 = np.concatenate([self.ux.flatten(),
                             first_qh_node[:,0].flatten(),
                             second_qh_node[:,0].flatten(),])
        vmin1 = np.min(mm1)
        vmax1 = np.max(mm1)

        mm2 = np.concatenate([self.uy.flatten(),
                             first_qh_node[:,1].flatten(),
                             second_qh_node[:,1].flatten()])

        vmin2 = np.min(mm2)
        vmax2 = np.max(mm2)

        import matplotlib.pyplot as plt
        # fig = plt.figure(dpi=300)
        # axs = fig.subplots(2,3, sharex=True, sharey=True, constrained_layout=True)
        # fig.subplots_adjust(hspace=0.1)
        levels = 10

        fig, axs = plt.subplots(nrows=2,ncols=3,dpi=300,sharex=True, sharey=True, figsize=(6,3))
        fig.subplots_adjust(wspace=0.1)
        tc1 = axs[0,0].tricontourf(self.x, self.y, self.triangles, self.ux, levels=levels, vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,0].tricontour(self.x, self.y, self.triangles, self.ux, linewidths=0.2, colors='k', levels=levels, vmin=vmin1, vmax=vmax1)

        axs[0,0].set_aspect('equal','box')
        axs[0,0].tick_params(bottom=False)
        # axs[0,0].tick_params(left=False, bottom=False)
        tc2 = axs[1,0].tricontourf(self.x, self.y, self.triangles, self.uy, levels=levels, vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,0].tricontour(self.x, self.y, self.triangles, self.uy, linewidths=0.2, colors='k', levels=levels, vmin=vmin2, vmax=vmax2)
        axs[1,0].set_aspect('equal','box')

        axs[0,1].tricontourf(first_x, first_y, first_tri, first_qh_node[:,0], levels=levels, vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,1].tricontour(first_x, first_y, first_tri, first_qh_node[:,0], linewidths=0.2, colors='k', levels=levels, vmin=vmin1, vmax=vmax1)
        axs[0,1].set_aspect('equal','box')
        axs[0,1].tick_params(left=False, bottom=False)

        axs[1,1].tricontourf(first_x, first_y, first_tri, first_qh_node[:,1], levels=levels, vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,1].tricontour(first_x, first_y, first_tri, first_qh_node[:,1], linewidths=0.2, colors='k', levels=levels, vmin=vmin2, vmax=vmax2)
        axs[1,1].set_aspect('equal','box')
        axs[1,1].tick_params(left=False)

        axs[0,2].tricontourf(second_x, second_y, second_tri, second_qh_node[:,0], levels=levels, vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,2].tricontour(second_x, second_y, second_tri, second_qh_node[:,0], linewidths=0.2, colors='k', levels=levels, vmin=vmin1, vmax=vmax1)
        axs[0,2].set_aspect('equal','box')
        axs[0,2].tick_params(left=False, bottom=False)

        axs[1,2].tricontourf(second_x, second_y, second_tri, second_qh_node[:,1], levels=levels, vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,2].tricontour(second_x, second_y, second_tri, second_qh_node[:,1], linewidths=0.2, colors='k', levels=levels, vmin=vmin2, vmax=vmax2)
        axs[1,2].set_aspect('equal','box')
        axs[1,2].tick_params(left=False)

        axs[0,0].set_yticks([0,0.5,1.0])
        axs[1,0].set_yticks([0,0.5,1.0])
        axs[1,0].set_xticks([0,0.5,1.0])
        axs[1,1].set_xticks([0,0.5,1.0])
        axs[1,2].set_xticks([0,0.5,1.0])


        fig.colorbar(tc1, ax=axs[0,:].flat)

        fig.colorbar(tc2, ax=axs[1,:].flat)

        # fig.tight_layout()
        plt.savefig(file_name)
        plt.show()
        # self.model.uh
        # self.model.qh




    def cmp_field_1(self, ux, uy, mesh_size, display=False, file_name='field.svg', **kwargs):
        self.model.Geometry = self.geometry_func(mesh_size)
        self.model.mesh(display=display)
        self.model.fem(display=display)
        self.model.efem(display=display)

        self.triangles = self.model.Element.copy()
        self.x = self.model.Node[:,0].copy()
        self.y = self.model.Node[:,1].copy()
        self.ux = ux(self.x, self.y)
        self.uy = uy(self.x, self.y)

        mm1 = np.concatenate([self.ux.flatten(),self.uy.flatten()])
        vmin1 = np.min(mm1)
        vmax1 = np.max(mm1)

        mm2 = np.concatenate([self.model.qh[:,0].flatten(),self.model.qh[:,1].flatten()])
        vmin2 = np.min(mm2)
        vmax2 = np.max(mm2)

        mm3 = np.concatenate([self.model.qh_fem[:,0].flatten(), self.model.qh_fem[:,1].flatten()])
        vmin3 = np.min(mm3)
        vmax3 = np.max(mm3)

        import matplotlib.pyplot as plt
        import matplotlib
        from matplotlib import ticker

        tick_locator = ticker.MaxNLocator(nbins=5)
        matplotlib.rc('text', usetex=True)
        matplotlib.rc('text.latex', preamble=r'\usepackage{amsmath}')

        fig = plt.figure(dpi=300)
        axs = fig.subplots(3,2, sharex=True, sharey=True)

        tc = axs[0,0].tricontourf(self.x, self.y, self.triangles, self.ux, levels=256, vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,0].set_aspect('equal','box')
        axs[0,0].tick_params(bottom=False)
        # axs[0,0].set_ylabel(r'$\boldsymbol{q}$', rotation=0)
        # axs[0,0].set_title(r'$x$')
        # axs[0,0].tick_params(left=False, bottom=False)
        axs[0,1].tricontourf(self.x, self.y, self.triangles, self.uy, levels=256, vmin=vmin1, vmax=vmax1, **kwargs)
        axs[0,1].set_aspect('equal','box')
        axs[0,1].tick_params(left=False, bottom=False)
        # axs[0,1].set_title(r'$y$')
        fc = fig.colorbar(tc, ax=axs[0,:].flat)
        # fc.ax.locator_params(nbins=3)
        # fc.set_ticks([0.4,0.6,0.8])
        # fc.ax.set_yticklabels(ticks)
        fc.update_ticks()


        tc = axs[1,0].tripcolor(self.x, self.y, self.triangles, facecolors=self.model.qh[:,0], vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,0].set_aspect('equal','box')
        # axs[1,0].set_ylabel(r'$\boldsymbol{q}_h$', rotation=0)
        axs[1,1].tripcolor(self.x, self.y, self.triangles, facecolors=self.model.qh[:,1], vmin=vmin2, vmax=vmax2, **kwargs)
        axs[1,1].set_aspect('equal','box')
        axs[1,1].tick_params(left=False)
        fc = fig.colorbar(tc, ax=axs[1,:].flat)
        fc.ax.locator_params(nbins=6)
        # fc.locator = tick_locator
        fc.update_ticks()


        tc = axs[2,0].tripcolor(self.x, self.y, self.triangles, facecolors=self.model.qh_fem[:,0], vmin=vmin3, vmax=vmax3, **kwargs)
        axs[2,0].set_aspect('equal','box')
        # axs[2,0].set_ylabel(r'$\nabla u_h$', rotation=0)
        axs[2,1].tripcolor(self.x, self.y, self.triangles, facecolors=self.model.qh_fem[:,1], vmin=vmin3, vmax=vmax3, **kwargs)
        axs[2,1].set_aspect('equal','box')
        axs[2,1].tick_params(left=False)
        fc = fig.colorbar(tc, ax=axs[2,:].flat)
        fc.ax.locator_params(nbins=5)
        fc.update_ticks()


        axs[0,0].set_yticks([0,0.5,1.0])
        axs[1,0].set_yticks([0,0.5,1.0])
        axs[2,0].set_yticks([0,0.5,1.0])

        # fig.tight_layout()
        plt.savefig(file_name)
        plt.show()
        # self.model.uh
        # self.model.qh

    def _cmp_point_over_mehes(self, start, end, num, pt, q):
        meshes = np.logspace(start, end, num)
        print("size\tqx(efem)\tqx(fem)\tqy(efem)\tqy(fem)\taeqx(efem)\taeqx(fem)\treqy(efem)\treqy(fem)")
        print("------------------------------------")

        for mesh_size in meshes:
            with MutePrints():
                self.model.Geometry = self.geometry_func(mesh_size)
                self.model.mesh()
                self.model.fem()
                self.model.efem()

            loc_node = np.where((self.model.Node[:,0]==pt[0])*(self.model.Node[:,1]==pt[1]))[0][0]
            # print(loc_node)
            if loc_node:
                loc_fem = []
                loc_efem = []

                for holeid, nodes in enumerate(self.model.Element):
                    if loc_node in nodes:
                        loc_efem.append(self.model.qh[holeid])
                        loc_fem.append(self.model.qh_fem[holeid])
                qh_at_pt = np.mean(np.array(loc_efem),axis=0)
                qh_fem_at_pt = np.mean(np.array(loc_fem),axis=0)

                print('{:.3f}\t{:.5f}\t{:.5f}\t{:.5f}\t{:.5f} \
                    \t{:.5f}\t{:.5f}\t{:.4%}\t{:.4%}'.format(mesh_size,
                        qh_at_pt[0],qh_fem_at_pt[0], qh_at_pt[1],qh_fem_at_pt[1],
                        qh_at_pt[0]-q[0], qh_fem_at_pt[0]-q[0], np.abs((qh_at_pt[1]-q[1])/q[1]),np.abs((qh_fem_at_pt[1]-q[1])/q[1])))

    def cmp_cre_field(self, ux, uy, mesh_size, display=False, file_name='field.svg', **kwargs):
        self.model.Geometry = self.geometry_func(mesh_size)
        self.model.mesh(display=display)
        self.model.fem(display=display)
        self.model.efem(display=display)

        self.triangles = self.model.Element.copy()
        self.x = self.model.Node[:,0].copy()
        self.y = self.model.Node[:,1].copy()
        self.ux = ux(self.x, self.y)
        self.uy = uy(self.x, self.y)

        mm = np.concatenate([(self.model.qh_fem[:,0]-self.model.qh[:,0]).flatten(), (self.model.qh_fem[:,1]-self.model.qh[:,1]).flatten()])
        vmin = np.min(mm)
        vmax = np.max(mm)

        import matplotlib.pyplot as plt
        fig = plt.figure(dpi=300)
        ax = fig.subplots(1,sharex=True, sharey=True)

        tc = ax.tripcolor(self.x, self.y, self.triangles, facecolors=self.model.qh_fem[:,0]-self.model.qh[:,0], vmin=vmin, vmax=vmax, **kwargs)
        ax.set_aspect('equal','box')

        ax.tripcolor(self.x, self.y, self.triangles, facecolors=self.model.qh_fem[:,1]-self.model.qh[:,1], vmin=vmin, vmax=vmax, **kwargs)
        ax.set_aspect('equal','box')
        ax.tick_params(left=False)

        ax.set_yticks([0,0.5,1.0])

        fig.colorbar(tc, ax=ax)

        # fig.tight_layout()
        plt.savefig(file_name)
        plt.show()
        # self.model.uh
        # self.model.qh

    def _cmp_field(self, u, ux, uy, mesh_size, display=False, file_name='field.svg'):
        self.model.Geometry = self.geometry_func(mesh_size)
        self.model.mesh(display=display)
        self.model.fem()
        self.model.efem()

        self.triangles = self.model.Element.copy()
        self.x = self.model.Node[:,0].copy()
        self.y = self.model.Node[:,1].copy()
        self.u = u(self.x, self.y)
        self.ux = ux(self.x, self.y)
        self.uy = uy(self.x, self.y)


        import matplotlib.pyplot as plt
        fig = plt.figure(dpi=300)
        axs = fig.subplots(2,3, sharex=True, sharey=True)
        def fmt(x):
            s = f"{x:.2f}"
            return rf"{s}" if plt.rcParams["text.usetex"] else f"{s}"

        tc = axs[0,0].tricontour(self.x, self.y, self.triangles, self.u, linewidths=0.5, colors='k', levels=15)
        tf = axs[0,0].tricontourf(self.x, self.y, self.triangles, self.u, levels=15)
        axs[0,0].clabel(tc, tc.levels, inline=True, fmt=fmt, fontsize=6, colors='white')
        axs[0,0].set_aspect('equal','box')
        axs[0,0].tick_params(bottom=False)

        tc = axs[1,0].tricontour(self.x, self.y, self.triangles, self.model.uh, linewidths=0.5, colors='k', levels=15)
        tf = axs[1,0].tricontourf(self.x, self.y, self.triangles, self.model.uh, levels=15)
        axs[1,0].clabel(tc, tc.levels, inline=True, fmt=fmt, fontsize=6, colors='white')
        axs[1,0].set_aspect('equal','box')

        tc = axs[0,1].tricontourf(self.x, self.y, self.triangles, self.ux, levels=256)
        axs[0,1].set_aspect('equal','box')
        axs[0,1].tick_params(left=False, bottom=False)

        axs[1,1].tripcolor(self.x, self.y, self.triangles, facecolors=self.model.qh[:,0])
        axs[1,1].set_aspect('equal','box')
        axs[1,1].tick_params(left=False)

        axs[0,2].tricontourf(self.x, self.y, self.triangles, self.uy, levels=256)
        axs[0,2].set_aspect('equal','box')
        axs[0,2].tick_params(left=False, bottom=False)

        axs[1,2].tripcolor(self.x, self.y, self.triangles, facecolors=self.model.qh[:,1])
        axs[1,2].set_aspect('equal','box')
        axs[1,2].tick_params(left=False)

        # plt.subplots_adjust(bottom=0.1, right=0.8, top=0.9)
        # newax = plt.axes([0.85, 0.1, 0.075, 0.9])
        fig.colorbar(tc, ax=axs.ravel().tolist())
        # tc = axs[0,0].tricontour(self.x, self.y, self.triangles, self.u, linewidths=0.5, colors='k', levels=15)
        # tf = axs[0,0].tricontourf(self.x, self.y, self.triangles, self.u, levels=15)
        # axs[0,0].clabel(tc, tc.levels, inline=True, fmt=fmt, fontsize=6, colors='white')
        # axs[0,0].set_aspect('equal','box')

        # tc = axs[1,0].tricontour(self.x, self.y, self.triangles, self.model.uh, linewidths=0.5, colors='k', levels=15)
        # tf = axs[1,0].tricontourf(self.x, self.y, self.triangles, self.model.uh, levels=15)
        # axs[1,0].clabel(tc, tc.levels, inline=True, fmt=fmt, fontsize=6, colors='white')
        # axs[1,0].set_aspect('equal','box')


        # tc = axs[0].tricontour(self.x, self.y, self.triangles, self.model.uh, linewidths=0.5, colors='k', levels=15)
        # tf = axs[0].tricontourf(self.x, self.y, self.triangles, self.model.uh, levels=15)
        # axs[0].clabel(tc, tc.levels, inline=True, fmt=fmt, fontsize=6)
        # axs[0].set_aspect('equal','box')


        # fig.colorbar(im, ax=axs[3])


        # fig.tight_layout()
        plt.savefig(file_name)
        plt.show()
        # self.model.uh
        # self.model.qh

    def __cmp_field(self, u, ux, uy, mesh_size, display=False):
        self.model.Geometry = self.geometry_func(mesh_size)
        self.model.mesh(display=display)
        self.model.fem()
        self.model.efem()

        self.triangles = self.model.Element.copy()
        self.x = self.model.Node[:,0].copy()
        self.y = self.model.Node[:,1].copy()


        import matplotlib.pyplot as plt
        fig = plt.figure(dpi=300)
        axs = fig.subplots(1,3)
        axs[0].triplot(self.x, self.y, self.triangles)
        axs[0].set_aspect('equal','box')
        axs[0].set_xlim([0,1])
        axs[0].set_ylim([0,1])


        axs[1].tripcolor(self.x, self.y, self.triangles, self.model.uh, shading='flat', edgecolors='black', linewidths=1)
        axs[1].set_aspect('equal','box')
        axs[1].set_xlim([0,1])
        axs[1].set_ylim([0,1])
        def fmt(x):
            s = f"{x:.2f}"
            # if s.endswith("0"):
            #     s = f"{x:.0f}"
            return rf"{s}" if plt.rcParams["text.usetex"] else f"{s}"

        tc = axs[2].tricontour(self.x, self.y, self.triangles, self.model.uh, linewidths=0.5, colors='k', levels=10)
        tf = axs[2].tricontourf(self.x, self.y, self.triangles, self.model.uh, levels=10)
        # print(type(tc), type(tf))
        axs[2].clabel(tc, tc.levels, inline=True, fmt=fmt, fontsize=6)
        axs[2].set_aspect('equal','box')
        # fig.colorbar(im, ax=axs[3])
        fig.tight_layout()
        plt.show()
        # self.model.uh
        # self.model.qh

    def cmp_flux(self, u, ux, uy, mesh_size, display=False, file_name='flux_field{}.svg'):
        self.model.Geometry = self.geometry_func(mesh_size)
        self.model.mesh(display=display)
        self.model.fem()
        self.model.efem()

        self.triangles = self.model.Element.copy()

        self.x = np.zeros(len(self.triangles))
        self.y = np.zeros(len(self.triangles))
        self.qx = np.zeros(len(self.triangles))
        self.qy = np.zeros(len(self.triangles))


        self.u = u(self.x, self.y)
        self.ux = ux(self.x, self.y)
        self.uy = uy(self.x, self.y)


        qh_fem = np.zeros((len(self.model.Element),2))

        for elemid, e in enumerate(self.model.Element):
            x1, y1 = self.model.Node[e][0]
            x2, y2 = self.model.Node[e][1]
            x3, y3 = self.model.Node[e][2]
            
            b1 = y2 - y3
            b2 = y3 - y1
            b3 = y1 - y2 
            c1 = x3 - x2
            c2 = x1 - x3 
            c3 = x2 - x1 

            DerPhiKe = np.array([[b1, b2, b3],
                                 [c1, c2, c3]])/(2*self.model.area[elemid])

            # construct the flux from EFEM solution
            qh_fem[elemid] = DerPhiKe @ self.model.uh[e]

            # construct the true flux field
            self.x[elemid], self.y[elemid] = np.mean(self.model.Node[e], axis=0)
            self.qx[elemid] = ux(self.x[elemid], self.y[elemid])
            self.qy[elemid] = uy(self.x[elemid], self.y[elemid])


        import matplotlib.pyplot as plt

        fig, axs = plt.subplots(1,3,dpi=300,sharey=True)
        # axs[0].triplot(self.model.Node[:,0], self.model.Node[:,1], self.triangles, linewidth=0.2, color='black', alpha=0.2)
        axs[0].quiver(self.x, self.y, self.qx, self.qy, angles='xy', color='blue',width=0.005, headwidth=3., headlength=4.)
        axs[0].set_aspect('equal', 'box')
        # axs[0].axis('off')
        
        # axs[1].triplot(self.model.Node[:,0], self.model.Node[:,1], self.triangles, linewidth=0.2, color='black', alpha=0.2)
        axs[1].quiver(self.x, self.y, qh_fem[:,0], qh_fem[:,1], angles='xy', color='blue',width=0.005, headwidth=3., headlength=4.)
        axs[1].set_aspect('equal', 'box')
        # axs[1].axis('off')

        # axs[2].triplot(self.model.Node[:,0], self.model.Node[:,1], self.triangles, linewidth=0.2, color='black', alpha=0.2)
        axs[2].quiver(self.x, self.y, self.model.qh[:,0], self.model.qh[:,1], angles='xy', color='blue',width=0.005, headwidth=3., headlength=4.)
        axs[2].set_aspect('equal', 'box')
        # axs[2].axis('off')
        fig.tight_layout()
        plt.savefig(file_name.format(mesh_size))
        plt.show()

    def cmp_mesh(self, mesh_size_list, display=False, file_name='Mesh{}.svg', color=None):
        
        import matplotlib.pyplot as plt
        for mesh_size in mesh_size_list:
            self.model.Geometry = self.geometry_func(mesh_size)
            self.model.mesh(display=display)
            print(self.model.numNode, 1/np.sqrt(self.model.numNode), mesh_size)
            self.triangles = self.model.Element.copy()
            self.x = self.model.Node[:,0].copy()
            self.y = self.model.Node[:,1].copy()

            fig, ax = plt.subplots(dpi=300)

            ax.triplot(self.x, self.y, self.triangles, color=color)
            ax.set_aspect('equal', 'box')
            ax.axis('off')
            fig.tight_layout()
            plt.savefig(file_name.format(mesh_size))

    def cmp_cond(self, cr, fem_kwargs={}, efem_kwargs={}, file_name='Cond{}.svg', font_size=18):

        A = np.c_[np.log(cr.hs_efem), np.ones_like(cr.hs_efem)]
        y1 = np.log(cr.conds_efem)
        y2 = np.log(cr.conds_fem)
        slope1, _ = np.linalg.inv(A.T@A) @ (A.T@y1)   
        slope2, _ = np.linalg.inv(A.T@A) @ (A.T@y2)   
        efem_kwargs['label'] = efem_kwargs['label'] + '(Slope: {:.2f})'.format(slope1)
        fem_kwargs['label'] = fem_kwargs['label'] + '(Slope: {:.2f})'.format(slope2)

        hs_group = [cr.hs_fem, cr.hs_efem]
        conds_group = [cr.conds_fem, cr.conds_efem]

        kwargs_group = [fem_kwargs, efem_kwargs] 

        import matplotlib
        import matplotlib.pyplot as plt
        matplotlib.rc('text', usetex=True)
        matplotlib.rc('text.latex', preamble=r'\usepackage{amsmath}')

        fig, ax = plt.subplots(dpi=300)
        ax.set_box_aspect(1)

        for h, e, kwargs in zip(hs_group, conds_group, kwargs_group):
            ax.plot(h, e, **kwargs) 

        ax.tick_params(axis="both", which="both", direction="in")    
        ax.tick_params(axis='both', which='major', labelsize=font_size)
        ax.tick_params(axis='both', which='major', length=8)
        ax.tick_params(axis='both', which='major', width=1)
        ax.tick_params(axis='both', which='minor', length=4)
        ax.tick_params(axis='both', which='minor', width=1)
        ax.set_xscale("log")
        ax.set_yscale("log")               
        ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())

        ax.xaxis.set_minor_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.LogFormatter(labelOnlyBase=False, minor_thresholds=(0.5,0.4)))
        plt.setp(ax.spines.values(), linewidth=2)
        plt.legend(fontsize='large', frameon=False)
        plt.tight_layout()
        plt.savefig(file_name)
        plt.show()

    def cmp_flux_error(self, ux, uy, mesh_size, display=False, file_name='flux_error.svg'):
        self.model.Geometry = self.geometry_func(mesh_size)
        self.model.mesh(display=display)
        self.model.fem()
        self.model.error_fem(ux, uy)
        self.model.efem()
        self.model.error_efem(ux, uy)

        self.triangles = self.model.Element.copy()
        self.x = self.model.Node[:,0].copy()
        self.y = self.model.Node[:,1].copy()


        norm_fem = np.sqrt(np.sum(self.model.flux_re_fem))
        norm_efem = np.sqrt(np.sum(self.model.flux_re_efem))


        import matplotlib.pyplot as plt
        fig = plt.figure(dpi=300)
        axs = fig.subplots(1,2, sharex=True, sharey=True)

        tc = axs[0].tripcolor(self.x, self.y, self.triangles, facecolors=np.sqrt(self.model.flux_re_fem)/norm_fem )
        axs[0].set_aspect('equal','box')
        axs[0].tick_params(left=False)

        axs[1].tripcolor(self.x, self.y, self.triangles, facecolors=np.sqrt(self.model.flux_re_efem)/norm_efem)
        axs[1].set_aspect('equal','box')
        axs[1].tick_params(left=False)

        # plt.subplots_adjust(bottom=0.1, right=0.8, top=0.9)
        # newax = plt.axes([0.85, 0.1, 0.075, 0.9])
        fig.colorbar(tc, ax=axs.ravel().tolist())
        # tc = axs[0,0].tricontour(self.x, self.y, self.triangles, self.u, linewidths=0.5, colors='k', levels=15)
        # tf = axs[0,0].tricontourf(self.x, self.y, self.triangles, self.u, levels=15)
        # axs[0,0].clabel(tc, tc.levels, inline=True, fmt=fmt, fontsize=6, colors='white')
        # axs[0,0].set_aspect('equal','box')

        # tc = axs[1,0].tricontour(self.x, self.y, self.triangles, self.model.uh, linewidths=0.5, colors='k', levels=15)
        # tf = axs[1,0].tricontourf(self.x, self.y, self.triangles, self.model.uh, levels=15)
        # axs[1,0].clabel(tc, tc.levels, inline=True, fmt=fmt, fontsize=6, colors='white')
        # axs[1,0].set_aspect('equal','box')


        # tc = axs[0].tricontour(self.x, self.y, self.triangles, self.model.uh, linewidths=0.5, colors='k', levels=15)
        # tf = axs[0].tricontourf(self.x, self.y, self.triangles, self.model.uh, levels=15)
        # axs[0].clabel(tc, tc.levels, inline=True, fmt=fmt, fontsize=6)
        # axs[0].set_aspect('equal','box')


        # fig.colorbar(im, ax=axs[3])

        # fig.tight_layout()
        plt.savefig(file_name)
        plt.show()
        # self.model.uh
        # self.model.qh

    def cmp_convergence_rate(self, cr_model, which='size', file_name=None, font_size=18, fem_kwargs={}, efem_kwargs={}, cre_kwargs={}):
        if (cr_model.hs_efem is None) or (cr_model.hs_fem is None):
            print("Please call convergence_all() method first!")
            return

        hs_group = [cr_model.hs_fem, cr_model.hs_efem]
        dofs_group = [cr_model.dofs_fem, cr_model.dofs_efem]
        errs_group = [cr_model.errs_fem, cr_model.errs_efem]
        dofs_slope_group = [cr_model.dof_slope_fem, cr_model.dof_slope_efem]
        size_slope_group = [cr_model.size_slope_fem, cr_model.size_slope_efem]

        kwargs_group = [fem_kwargs, efem_kwargs] 

        import matplotlib.pyplot as plt
        import matplotlib

        matplotlib.rc('text', usetex=True)
        matplotlib.rc('text.latex', preamble=r'\usepackage{amsmath}')

        fig, ax = plt.subplots(dpi=300)
        # ax.set_aspect('equal', 'box')
        ax.set_box_aspect(1)
        if which == "size":
            for h, e, s, kwargs in zip(hs_group, errs_group, size_slope_group, kwargs_group):
                kwargs['label'] = kwargs['label'] + '(Slope: {:.2f})'.format(s)
                ax.plot(h, e, **kwargs)

            cre_kwargs['label'] = cre_kwargs['label'] + '(Slope: {:.2f})'.format(cr_model.size_slope_cre)
            ax.plot(cr_model.hs_fem, cr_model.cres, **cre_kwargs)

        elif which == "dofs":
            for d, e, s, kwargs in zip(dofs_group, errs_group, dofs_slope_group, kwargs_group):
                kwargs['label'] = kwargs['label'] + '(Slope: {:.2f})'.format(s)
                ax.plot(d, e, **kwargs)
            
            cre_kwargs['label'] = cre_kwargs['label'] + '(Slope: {:.2f})'.format(cr_model.dof_slope_cre)
            ax.plot(cr_model.dofs_fem, cr_model.cres, **cre_kwargs)
        else:
            print("The parameter 'which' should be 'size' or 'dofs'!")
            return    

        ax.tick_params(axis="both", which="both", direction="in")    
        ax.tick_params(axis='both', which='major', labelsize=font_size)
        ax.tick_params(axis='both', which='major', length=8)
        ax.tick_params(axis='both', which='major', width=1)
        ax.tick_params(axis='both', which='minor', length=4)
        ax.tick_params(axis='both', which='minor', width=1)

        # ax.tick_params(axis = 'both', which = 'minor', labelsize = 24)

        ax.set_xscale("log")
        ax.set_yscale("log")
        # ax.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        
        # ax.set_ylabel(r"Error in norm, $e$")
        # ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        # ax.xaxis.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%d'))
        ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.LogFormatter(labelOnlyBase=True, minor_thresholds=(0.5,0.2)))
        
        # ax.set_xlim([0.0085,0.04])
        # ax.set_ylim([0.07, 1.1])
        # ax.set_xlim([0.005,0.5])
        # ax.set_ylim([0.05,20])
        # ax.grid()
        plt.setp(ax.spines.values(), linewidth=2)
        plt.legend(fontsize='large',frameon=False)
        plt.tight_layout()
        if file_name is None:
            plt.savefig("multiple_convergence.png")
        else:
            plt.savefig(file_name)
        plt.show()        

    def cre_convergence_rate(self, cr_model, which='size', file_name=None, font_size=18, cre_kwargs={}):
        if (cr_model.hs_efem is None) or (cr_model.hs_fem is None):
            print("Please call convergence_all() method first!")
            return

        import matplotlib.pyplot as plt
        import matplotlib

        matplotlib.rc('text', usetex=True)
        matplotlib.rc('text.latex', preamble=r'\usepackage{amsmath}')

        fig, ax = plt.subplots(dpi=300)
        ax.set_box_aspect(1)
        # ax.set_aspect('equal', 'box')
        if which == "size":
            cre_kwargs['label'] = cre_kwargs['label'] + '(Slope: {:.2f})'.format(cr_model.size_slope_cre)
            ax.plot(cr_model.hs_fem, cr_model.cres, **cre_kwargs)

        elif which == "dofs":
            cre_kwargs['label'] = cre_kwargs['label'] + '(Slope: {:.2f})'.format(cr_model.dof_slope_cre)
            ax.plot(cr_model.dofs_fem, cr_model.cres, **cre_kwargs)
        else:
            print("The parameter 'which' should be 'size' or 'dofs'!")
            return    

        ax.tick_params(axis="both", which="both", direction="in")    
        ax.tick_params(axis='both', which='major', labelsize=font_size)
        ax.tick_params(axis='both', which='major', length=8)
        ax.tick_params(axis='both', which='major', width=1)
        ax.tick_params(axis='both', which='minor', length=4)
        ax.tick_params(axis='both', which='minor', width=1)

        # ax.tick_params(axis = 'both', which = 'minor', labelsize = 24)

        ax.set_xscale("log")
        ax.set_yscale("log")
        # ax.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        
        # ax.set_ylabel(r"Error in norm, $e$")
        # ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        # ax.xaxis.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%d'))
        ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.LogFormatter(labelOnlyBase=True, minor_thresholds=(0.5,0.2)))
        ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
        # ax.set_ylim([0.2,2])
        # ax.grid()
        plt.setp(ax.spines.values(), linewidth=2)
        plt.legend(fontsize='large', frameon=False)
        if file_name is None:
            plt.savefig("multiple_convergence.png")
        else:
            plt.savefig(file_name)
        plt.show()        

    def resultant_force_convergence_rate(self, cr_model, q0n, file_name="r1err.png", font_size=18, kwargs=[]):

        import matplotlib.pyplot as plt
        import matplotlib

        matplotlib.rc('text', usetex=True)
        matplotlib.rc('text.latex', preamble=r'\usepackage{amsmath}')

        fig, ax = plt.subplots(dpi=300)
        ax.set_box_aspect(1)

        for holeid in range(cr_model.model.numHole):

            rs = np.array([x[0] for x in cr_model.rs])
            rs_fem = np.array(cr_model.rs_fem)
            err = np.abs(rs[:,holeid]-rs_fem[:,holeid]+q0n[holeid])

            h = 1/np.sqrt(cr_model.dofs_fem)

            A = np.c_[np.log(h), np.ones_like(h)]
            y = np.log(err)
            k, _ = np.linalg.inv(A.T@A) @ (A.T@y)

            kwargs[holeid]['label'] = kwargs[holeid]['label'] + '(Slope: {:.2f})'.format(k)
            ax.plot(h, err, **kwargs[holeid])


        ax.tick_params(axis="both", which="both", direction="in")    
        ax.tick_params(axis='both', which='major', labelsize=font_size)
        ax.tick_params(axis='both', which='major', length=8)
        ax.tick_params(axis='both', which='major', width=1)
        ax.tick_params(axis='both', which='minor', length=4)
        ax.tick_params(axis='both', which='minor', width=1)

        # ax.tick_params(axis = 'both', which = 'minor', labelsize = 24)

        ax.set_xscale("log")
        ax.set_yscale("log")
        # ax.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        
        # ax.set_ylabel(r"Error in norm, $e$")
        # ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        # ax.xaxis.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%d'))
        ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.LogFormatter(labelOnlyBase=True, minor_thresholds=(0.5,0.2)))
        ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
        # ax.set_ylim([0.2,2])
        # ax.grid()
        plt.setp(ax.spines.values(), linewidth=2)
        plt.legend(fontsize='large', frameon=False)

        plt.savefig(file_name)
        plt.show()        

    def _resultant_force_convergence_rate(self, cr_model, q0n, holeid=0, file_name="r1err.png", font_size=18, **kwargs):

        import matplotlib.pyplot as plt
        import matplotlib

        matplotlib.rc('text', usetex=True)
        matplotlib.rc('text.latex', preamble=r'\usepackage{amsmath}')

        fig, ax = plt.subplots(dpi=300)
        ax.set_box_aspect(1)


        
        rs = np.array([x[0] for x in cr_model.rs])
        rs_fem = np.array(cr_model.rs_fem)
        err = np.abs(rs[:,holeid]-rs_fem[:,holeid]+q0n[holeid])

        h = 1/np.sqrt(cr_model.dofs_fem)

        A = np.c_[np.log(h), np.ones_like(h)]
        y = np.log(err)
        k, _ = np.linalg.inv(A.T@A) @ (A.T@y)


        kwargs['label'] = kwargs['label'] + '(Slope: {:.2f})'.format(k)
        ax.plot(h, err, **kwargs)


        ax.tick_params(axis="both", which="both", direction="in")    
        ax.tick_params(axis='both', which='major', labelsize=font_size)
        ax.tick_params(axis='both', which='major', length=8)
        ax.tick_params(axis='both', which='major', width=1)
        ax.tick_params(axis='both', which='minor', length=4)
        ax.tick_params(axis='both', which='minor', width=1)

        # ax.tick_params(axis = 'both', which = 'minor', labelsize = 24)

        ax.set_xscale("log")
        ax.set_yscale("log")
        # ax.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
        
        # ax.set_ylabel(r"Error in norm, $e$")
        # ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        # ax.xaxis.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%d'))
        ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.ScalarFormatter())
        ax.xaxis.set_minor_formatter(matplotlib.ticker.LogFormatter(labelOnlyBase=True, minor_thresholds=(0.5,0.2)))
        ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
        # ax.set_ylim([0.2,2])
        # ax.grid()
        plt.setp(ax.spines.values(), linewidth=2)
        plt.legend(fontsize='large', frameon=False)

        plt.savefig(file_name)
        plt.show()        




def show_cross_element(model):
    import numpy as np
    from mayavi import mlab

    x = model.Node[:,0]
    y = model.Node[:,1]
    z = np.zeros_like(x)
    triangles1 = model.Element
    triangles2 = model.Element[list(model.crossElement.keys()),:]

    mlab.triangular_mesh(x, y, z, triangles1, color=(1,1,1),representation='wireframe' )
    mlab.triangular_mesh(x, y, z, triangles2, color=(1,0,0),representation='surface' )

    for cluster in model.NeuCluster:
        i,_ = cluster[0]
        x = model.Node[i][0]
        y = model.Node[i][1]
        z = 0    
        mlab.points3d(x,y,z,1, color=(0,1,0),scale_factor=0.01)
        for i, j in cluster:
            x = model.Node[i][0]
            y = model.Node[i][1]
            z = 0
            u = model.Node[j][0] - x
            v = model.Node[j][1] - y
            w = 0
            mlab.quiver3d(x, y, z, u, v, w, scale_factor=1)

    for edge in model.crossEdge:
        i,j = edge
        x = [model.Node[i][0], model.Node[j][0]]
        y = [model.Node[i][1], model.Node[j][1]]
        z = [0,0]
        mlab.plot3d(x, y, z, color=(0,0,1),tube_radius=0.005)
            
            
    mlab.view(0,0)
    # print(mlab.view())
    mlab.show()


def display_flux(model):
    import numpy as np
    from mayavi import mlab

    x = model.Node[:,0]
    y = model.Node[:,1]
    z = np.zeros_like(x)
    triangles1 = model.Element

    mlab.triangular_mesh(x, y, z, triangles1, color=(1,1,1),representation='surface' )
    # mlab.triangular_mesh(x, y, z, triangles2, color=(1,0,0),representation='surface' )

    # for cluster in model.NeuCluster:
    #     i,_ = cluster[0]
    #     x = model.Node[i][0]
    #     y = model.Node[i][1]
    #     z = 0    
    #     mlab.points3d(x,y,z,1, color=(0,1,0),scale_factor=0.01)
    #     for i, j in cluster:
    #         x = model.Node[i][0]
    #         y = model.Node[i][1]
    #         z = 0
    #         u = model.Node[j][0] - x
    #         v = model.Node[j][1] - y
    #         w = 0
    #         mlab.quiver3d(x, y, z, u, v, w, scale_factor=1)

    # for edge in model.crossEdge:
    #     i,j = edge
    #     x = [model.Node[i][0], model.Node[j][0]]
    #     y = [model.Node[i][1], model.Node[j][1]]
    #     z = [0,0]
    #     mlab.plot3d(x, y, z, color=(0,0,1),tube_radius=0.005)
            
            
    mlab.view(0,0)
    # print(mlab.view())
    mlab.show()